/*
 * Grid.cpp
 *
 *  Created on: Jun 6, 2011
 *      Author: marchi
 */

#include "Grid.h"
matrix Grid::co={{0,0,0},{0,0,0},{0,0,0}};
matrix Grid::oc={{0,0,0},{0,0,0},{0,0,0}};
double Grid::Volume=0.0;
unsigned int Grid::nnx=0;
unsigned int Grid::nny=0;
unsigned int Grid::nnz=0;
Grid::Grid(){
		try{
			if(nnx+nny+nnz==0) throw "Should initialize Metrics first";
		}
		catch(const char * s){
			std::cout << s << std::endl;
			exit(1);
		}
}
Grid::Grid(Grid & y){
	Allocate();
	*this=y;
}
Grid::Grid(const Grid & y){
	Allocate();
	*this=y;
}

Grid::Grid(array3<Complex> & x){
	this->Allocate();
	*this=x;
}
Grid::Grid(array3<double> & x){
	this->Allocate();
	*this=x;
}
Grid::Grid(const double & x){
	this->Allocate();
	*this=x;
}

Grid & Grid::operator=(const double & y){
	double * pt=&(*this)[0][0][0];
	for(int i=0;i<static_cast<int>(nnx*nny*nnz);i++) pt[i]=y;
	return *this;
}
Grid & Grid::operator=(array3<double> & a){
	try{
		if(this->Size() != a.Size()) throw "Grid and array3 arrays do not have the same size ";
	}
	catch(char * s){
		std::cout << s << std::endl;
		exit(1);
	}
	for(unsigned int i=0;i<nnx;i++)
		(*this)[i]=a[i];
	return *this;
}
Grid & Grid::operator=(array3<Complex> & a){
	try{
		if(this->Size() != a.Size()) throw "Grid and array3 arrays do not have the same size ";
	}
	catch(char * s){
		std::cout << s << std::endl;
		exit(1);
	}
	for(unsigned i=0;i<nnx;i++)
		for(unsigned j=0;j<nny;j++)
			for(unsigned k=0;k<nnz;k++){
				(*this)[i][j][k]=a[i][j][k].real();
			}
	return *this;
}
Grid Grid::operator-(Grid & a){
	Grid temp=*this;
	temp-=a;
	return temp;
}
Grid Grid::operator+(Grid & a){
	Grid temp=*this;
	temp+=a;
	return temp;
}
Grid Grid::operator*(const double & a){
	Grid temp=*this;
	temp*=a;
	return temp;
}
Grid Grid::operator/(const double & a){
	Grid temp=*this;
	temp/=a;
	return temp;
}

void Grid::setMetric(const matrix & co_in)
{
  /* Save some time by assuming lower right part is zero */
  for(int i=0;i<DIM;i++)
	  for(int j=0;j<DIM;j++) co[i][j]=co_in[i][j];
  real tmp=1.0/(co[XX][XX]*co[YY][YY]*co[ZZ][ZZ]);
  oc[XX][XX]=co[YY][YY]*co[ZZ][ZZ]*tmp;
  oc[YY][XX]=0;
  oc[ZZ][XX]=0;
  oc[XX][YY]=-co[XX][YY]*co[ZZ][ZZ]*tmp;
  oc[YY][YY]=co[XX][XX]*co[ZZ][ZZ]*tmp;
  oc[ZZ][YY]=0;
  oc[XX][ZZ]=(co[XX][YY]*co[YY][ZZ]-co[YY][YY]*co[XX][ZZ])*tmp;
  oc[YY][ZZ]=-co[YY][ZZ]*co[XX][XX]*tmp;
  oc[ZZ][ZZ]=co[XX][XX]*co[YY][YY]*tmp;
  Volume=(co[XX][XX]*(co[YY][YY]*co[ZZ][ZZ]-co[ZZ][YY]*co[YY][ZZ])
		  -co[YY][XX]*(co[XX][YY]*co[ZZ][ZZ]-co[ZZ][YY]*co[XX][ZZ])
		  +co[ZZ][XX]*(co[XX][YY]*co[YY][ZZ]-co[YY][YY]*co[XX][ZZ]));
}

void Grid::Xplor(FILE * fp){
	int nx_s,ny_s,nz_s;
	int nx_e,ny_e,nz_e;

	double avg=0.0, avg2=0.0;
	double a=0.0,b=0.0,c=0.0,alpha=0.0, beta=0, gamma=0;
	Grid & MyGrid=*this;

	int nx0=static_cast<int>(nnx);
	int ny0=static_cast<int>(nny);
	int nz0=static_cast<int>(nnz);
	for(int i=0; i< nx0; i++){
		for(int j=0; j< ny0; j++){
			for(int k=0; k< nz0; k++){

				avg+=MyGrid[i][j][k];
				avg2+=MyGrid[i][j][k]*MyGrid[i][j][k];
			}

		}
	}
	avg/=static_cast<double>(nx0*ny0*nz0);
	avg2/=static_cast<double>(nx0*ny0*nz0);

	nx_s=0;
	ny_s=0;
	nz_s=0;
	nx_e=nx0-1+nx_s;
	ny_e=ny0-1+ny_s;
	nz_e=nz0-1+nz_s;


	alpha=Alpha()*RAD2DEG;
	beta=Beta()*RAD2DEG;
	gamma=Gamma()*RAD2DEG;

	a=norm(co[XX]);
	b=norm(co[YY]);
	c=norm(co[ZZ]);
	fprintf(fp,"\n       2 !NTITLE\n") ;
	fprintf(fp," REMARKS Electrostatic Potential from GROMACS\n") ;
	fprintf(fp," REMARKS DATE: 2010-04-21\n") ;
	fprintf(fp,"%8d%8d%8d%8d%8d%8d%8d%8d%8d\n",nx,nx_s,nx_e,ny,ny_s,ny_e,nz,nz_s,nz_e);
	fprintf(fp,"%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",10*a,10*b,10*c,alpha,beta,gamma);
	fprintf(fp, "ZYX\n") ;
	int ja;
	for(int k=0;k<nz0;k++){
		fprintf(fp,"%8d\n",k);
		ja=0;
		for(int j=0;j<ny0;j++){
			for(int i=0;i<nx0;i++){
				ja++;
				fprintf(fp,"%12.5E",MyGrid[i][j][k]);
				if(!(ja%6)) fprintf(fp,"\n");
			}
		}
		fprintf(fp,"\n");
	}
	fprintf(fp,"%8d\n",-9999);

	double stddev=sqrt(avg2-avg*avg);
	fprintf(fp,"%12.4E %12.4E\n",avg,stddev);

}

const double Grid::Alpha(){
	double alpha;
	if (norm2(co[YY])*norm2(co[ZZ])!=0)
		alpha = acos(cos_angle(co[YY],co[ZZ]));
	else
		alpha = 0.5*M_PI;

	return alpha;
}
const double Grid::Beta(){
	double beta;
	if (norm2(co[XX])*norm2(co[ZZ])!=0)
		beta  = acos(cos_angle(co[XX],co[ZZ]));
	else
		beta  = 0.5*M_PI;
	return beta;
}
const double Grid::Gamma(){
	double gamma;
	if (norm2(co[XX])*norm2(co[YY])!=0)
		gamma = acos(cos_angle(co[XX],co[YY]));
	else
		gamma = 0.5*M_PI;
	return gamma;
}
void Grid::Rdf(FILE * fp,const double x[DIM],const double & cut, const double & dx){
	double xa[DIM];
	int nrdf=(int) (cut/dx)+1;

	double * rdf=new double [nrdf];
	int * irdf=new int [nrdf];
	for(int i=0;i<nrdf;i++) {
		rdf[i]=0.0;
		irdf[i]=0;
	}
	xa[XX] = oc[XX][XX] * x[XX] + oc[XX][YY] * x[YY] + oc[XX][ZZ] * x[ZZ];
	xa[YY] = oc[YY][XX] * x[XX] + oc[YY][YY] * x[YY] + oc[YY][ZZ] * x[ZZ];
	xa[ZZ] = oc[ZZ][XX] * x[XX] + oc[ZZ][YY] * x[YY] + oc[ZZ][ZZ] * x[ZZ];
	int nx0=(int) nnx;
	int ny0=(int) nny;
	int nz0=(int) nnz;
	double rx,ry,rz;
	double rnx=1.0/ (double)(nx0);
	double rny=1.0/ (double)(ny0);
	double rnz=1.0/ (double)(nz0);
	double xd[DIM],xc[DIM],dist;
	for(int i=0;i<nx0;i++){
		rx=(double) (i)*rnx;
		for(int j=0;j<ny0;j++){
			ry=(double) (j)*rny;
			for(int k=0;k<nz0;k++){
				rz=(double) (k)*rnz;
				xd[XX]=rx-xa[XX];
				xd[YY]=ry-xa[YY];
				xd[ZZ]=rz-xa[ZZ];
				xd[XX]=xd[XX]-rint(xd[XX]);
				xd[YY]=xd[YY]-rint(xd[YY]);
				xd[ZZ]=xd[ZZ]-rint(xd[ZZ]);

				xc[XX]=co[XX][XX]*xd[XX]+co[XX][YY]*xd[YY]+co[XX][ZZ]*xd[ZZ];
				xc[YY]=co[YY][XX]*xd[XX]+co[YY][YY]*xd[YY]+co[YY][ZZ]*xd[ZZ];
				xc[ZZ]=co[ZZ][XX]*xd[XX]+co[ZZ][YY]*xd[YY]+co[ZZ][ZZ]*xd[ZZ];
				dist=sqrt(xc[XX]*xc[XX]+xc[YY]*xc[YY]+xc[ZZ]*xc[ZZ]);
				if(dist < cut){
					int h=(int) (dist/dx);
					rdf[h]+=(*this)[i][j][k];
					irdf[h]++;
				}

			}
		}
	}
	for(int i=1;i<nrdf;i++){
		double bx=(double) (i)*dx;
		if(irdf[i]) {
			double fx=rdf[i]/(double) (irdf[i]);
			fprintf(fp," %f  %e \n",bx,fx);
		}
	}
	fprintf(fp,"& \n");
}


Grid::~Grid() {
	// TODO Auto-generated destructor stub
}
